
# Load collapse package
library(collapse)
# Add day-of-month
sub_dt[, `:=`(
  dom = begin_date %>% day()
)]
# Within each city-month-utility, add 
# 1. the first observed price
# 2. the median price
setkey(sub_dt, utility, city_ym, begin_date)
sub_dt[, `:=`(
  price_mrg_lag2_first = ffirst(price_mrg_lag2),
  spot_price_1week_lag2_first = ffirst(spot_price_1week_lag2),
  price_mrg_lag2_median = fmedian(price_mrg_lag2),
  spot_price_1week_lag2_median = fmedian(spot_price_1week_lag2)
), by = .(utility, city_ym)]

# Standard specification
# Marginal price
est_stnd_mrg = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)

# Shut down the border discontinuity
# Strategy: Interact 'utility' with city-YM FEs
# Marginal price
est_nodisc_mrg = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym^utility |
    log(price_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)

# Shut down the within-month-utility price variation
# Strategy: Interact Day-of-month with city-YM FEs
# Marginal price
est_dom_mrg = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym^dom |
    log(price_mrg_lag2) ~ spot_price_1week_lag2 + spot_price_1week_lag2:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)

# Shut down the within-month-utility price variation
# Strategy: Use the first (or median) price observed for each utility-city-YM
# Marginal price
est_first_mrg = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_mrg_lag2_first) ~ spot_price_1week_lag2_first + spot_price_1week_lag2_first:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)
est_median_mrg = feols(
  log(thm_day) ~
    I(bill_hdd/days) | 
    hh_id + city_ym |
    log(price_mrg_lag2_median) ~ spot_price_1week_lag2_median + spot_price_1week_lag2_median:utility,
  cluster = ~hh_id + price_cluster,
  data = sub_dt,
  mem.clean = TRUE,
  lean = FALSE
)

# Make a table with the results
list(
  est_stnd_mrg,
  est_nodisc_mrg,
  est_dom_mrg,
  est_first_mrg,
  est_median_mrg
) %>% etable()